for (mod in c("lasso","rf","nn","gbm")) {
  load(paste(modeldir,"/",
           country,
           "_",mod,"_",
           v,
           "_",
           rhs.group,
           "_mse",
           ".RData",
           sep = ""))
}
i = 0
ebma.results <- list()
lasso.fit <- list()
rf.fit <- list()
svm.fit <- list()
ada.fit <- list()
nn.fit <- list()
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Fit LASSO
  ebma.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- length(dta.past[,dv])
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("glmnet",
                                       "randomForest",
                                       "gbm",
                                       "nnet",
                                       "pROC")) %dopar% {   # Loop over CV runs
           set.seed(52*cv) 
           shuffle <- sample(1:N,
                             N,
                             replace=FALSE)
           pred.prob <- cbind(lasso.prediction = matrix(NA,nrow=N,ncol=1),
                              rf.prediction = matrix(NA,nrow=N,ncol=1),
                              gbm.prediction = matrix(NA,nrow=N,ncol=1),
                              nn.prediction = matrix(NA,nrow=N,ncol=1),
                              ebma.prediction = matrix(NA,nrow=N,ncol=1))
           colnames(pred.prob) <- c("lasso.prediction" ,
                                    "rf.prediction",
                                    "gbm.prediction" ,
                                    "nn.prediction" ,
                                    "ebma.prediction")
           sampleSplit = c(0,round((N/5*c(1:5))))
           for (f in 1:5) {
             # Remove constant variables from predictors
             train.RHS <- dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                        rhs]
             test.RHS <- dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                       rhs]
             not.constant <- apply(train.RHS,2,sd)!=0
             train.RHS <- train.RHS[,not.constant]
             test.RHS <- test.RHS[,not.constant]
             # Fit Lasso
             cv.fit<- glmnet(x = as.matrix(train.RHS),
                    y = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                dv]),
                    alpha = lasso_params$alpha,
                    family = lasso_params$family,
                    lambda = seq(10.9,1,-.1)^5*lasso.results[[i]]$lambda
             )
             opt.lam <- 100
             while (is.na(cv.fit$lambda[opt.lam])) {
               opt.lam = opt.lam - 1
             }
             pred.prob[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                       "lasso.prediction"]<- predict(cv.fit,
                                                   newx = as.matrix(test.RHS),
                                                   type = "response")[,opt.lam]
             # Fit GBM
             cv.fit = gbm.fit(x = train.RHS,
                              y = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                          dv]),
                              n.trees = gbm.results[[i]]$trees,
                              distribution = gbm_params$distribution,
                              interaction.depth = gbm_params$interaction.depth,
                              n.minobsinnode =gbm_params$n.minobsinnode,
                              verbose = gbm_params$verbose,
                              shrinkage = gbm_params$shrinkage,
                              bag.fraction=gbm_params$bag.fraction,
                              train.fraction=gbm_params$train.fraction)
             pred.prob[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                       "gbm.prediction"] <- predict(cv.fit,
                                    newdata = test.RHS,
                                    type = "response",
                                    n.trees = gbm.results[[i]]$trees)
             # Fit Random Forest
             cv.fit = randomForest(x = as.matrix(train.RHS),
                                   y = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                               dv]),
                                   ntree = rf_params$ntree,
                                   maxnodes = rf_params$maxnodes,
                                   type = rf_params$type,
                                   importance = rf_params$importance
             )
             pred.prob[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                       "rf.prediction"] <- predict(cv.fit,
                                                   newdata = test.RHS)
             # Find PC weights
             pcs.fit <- prcomp(train.RHS,center = TRUE, scale = TRUE)
             train.pcs <- predict(pcs.fit,
                                  train.RHS)[,1:min(nn_params$pcNum[[country]],
                                                    ncol(train.RHS))]
             test.pcs <- predict(pcs.fit,
                                 test.RHS)[,1:min(nn_params$pcNum[[country]],
                                                  ncol(train.RHS))]
             # Fit Neural Network
             cv.fit = nnet(x = train.pcs,
                           y = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                       dv]),
                           size = nn.results[[i]]$size, 
                           decay = nn_params$decay,
                           trace = nn_params$trace
             )
             pred.prob[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                       "nn.prediction"] <- predict(cv.fit,
                                     newdata = as.matrix(test.pcs))
           }
           weight.denom <- colSums(pred.prob[,c("lasso.prediction",
                                                "rf.prediction",
                                                "gbm.prediction",
                                                "nn.prediction")] 
                                   * dta[dta[,t.var] < year,
                                              dv] 
                                   + (1 - pred.prob[,c("lasso.prediction",
                                                       "rf.prediction",
                                                       "gbm.prediction",
                                                       "nn.prediction")]) 
                                   * (1 - dta[dta[,t.var] < year,
                                                   dv]))
           weight <- weight.denom/sum(weight.denom)
           pred.prob[,"ebma.prediction"] <- as.matrix(pred.prob[,c("lasso.prediction",
                                                                   "rf.prediction",
                                                                   "gbm.prediction",
                                                                   "nn.prediction")]) %*% as.matrix(weight)
           opt.dt <- mostaccDT(y=dta[dta[,t.var] < year,
                                          dv],
                               yhat=pred.prob[,"ebma.prediction"],
                               J=1000)
           opt.dtsens <- bestDT(y=dta[dta[,t.var] < year,
                                          dv],
                               yhat=pred.prob[,"ebma.prediction"],
                               J=1000)
           cv_results <- c(cv,weight,opt.dt,
                           opt.dtsens)
           names(cv_results) <- c("CV Run",
                                  "lasso.weight",
                                  "rf.weight",
                                  "gbm.weight",
                                  "nn.weight",
                                  "dt",
                                  "dtsens") 
           return(cv_results)
         }
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  set.seed(i)
  ebma.results[[i]] <- list()
  if (cv.runs>1) {
    ebma.results[[i]]$optim.dt <- mean(cv_results[,"dt"])   # Find optimal lambda
    ebma.results[[i]]$optim.dtsens <- mean(cv_results[,"dtsens"])
    ebma.results[[i]]$optim.wts <- colMeans(cv_results[,c("lasso.weight",
                                                      "rf.weight",
                                                      "gbm.weight",
                                                      "nn.weight")])
  } else {
    ebma.results[[i]]$optim.dt <- cv_results["dt"]
    ebma.results[[i]]$optim.dtsens <- cv_results["dtsens"]
    ebma.results[[i]]$optim.wts <- cv_results[c("lasso.weight",
                                                "rf.weight",
                                                "gbm.weight",
                                                "nn.weight")]
  }
  ebma.results[[i]]$oos.component.prob <- cbind("lasso" = lasso.results[[i]]$fit.oos,
                                                 "rf" = rf.results[[i]]$fit.oos,
                                                 "gbm" = gbm.results[[i]]$fit.oos,
                                                 "nn" = nn.results[[i]]$fit.oos)

  ebma.results[[i]]$fit.oos <- ebma.results[[i]]$oos.component.prob %*% as.matrix(ebma.results[[i]]$optim.wts)
  ebma.results[[i]]$actual.oos <- as.matrix(dta[dta[,t.var] == year,
                                                    dv])
  print(year)
}
save(ebma.results,file=paste(modeldir,"/",
                           country,
                           "_ebma_",
                           v,
                           "_",
                           rhs.group,
                           "_mse",
                           ".RData",
                           sep = ""))